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We investigate the thermodynamic and critical properties of an interact- 
. ing domain wall model which is derived from the triangular lattice antiferro- 
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^SJ . magnetic Ising model with the anisotropic nearest and next nearest neighbor 

, interactions. The model is equivalent to the general five-vertex model. Diag- 
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onalizing the transfer matrix exactly by the Bethe Ansatz method, we obtain 
the phase diagram displaying the commensurate and incommensurate (IC) 



^ , phases separated by the Pokrovsky-Talapov transitions. The phase diagram 

^ exhibits commensurate phases where the domain wall density q is locked at 

^ , the values of 0, 1/2 and 1. The IC phase is a critical state described by the 

H : 

Gaussian fixed point. The effective Gaussian coupling constant is obtained 
analytically and numerically for the IC phase using the finite size scaling 
predictions of the conformal field theory. It takes the value 1/2 in the non- 
interacting limit and also at the boundaries of g = or 1 phase and the value 
2 at the boundary of q = 1/2 phase, while it varies smoothly throughout the 
IC region. 
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I. INTRODUCTION 



There has been much interest in the two dimensional statistical mechanical systems 
which exhibit modulated phases on the periodic substrate [|l|. Among those systems are 
monolayers of physisorbed gas on solid surface which display incommensurate (IC) and 
commensurate (C) phases. In the domain wall description of IC phases domain walls 
separating commensurate patches are considered as the basic fluctuating degrees of freedom. 
The domain walls can be arranged either parallel to each other (striped domain wall) or 
in hexagonal pattern (honeycomb domain wall) depending on the domain wall crossing 
energy [^. The simplest type of commensurate-incommensurate (C-IC) transition is the 
Pokrovsky-Talapov (PT) transition 0] which describes the transition into striped IC phase. 
Here, the fluctuations of the striped domain wall cause an effective repulsive interaction 
between walls. The interaction varies as if / is the average distance between walls. Due 
to this repulsive interaction between domain walls, the C-IC transition to the striped IC 
phase is a continuous transition with the specific heat exponent a = 1/2 and the domain 
wall density displays a square root dependence on the chemical potential of domain wall if 
we approach the phase boundary from the incommensurate side. The theory is explicitly 
realized in fermion models of striped IC phases where domain walls are represented as world 
lines of fermions living in one dimensional chain. Free fermion model is also obtained as low 
temperature approximation to the ANNNI model [^]. In these models, the IC phase is a 
critical phase where the correlation functions decay by the power laws of the distance rather 
than by the exponential function of the distance. Recently, Park and Widom showed that 
the IC phase modeled by free fermion hamiltonian is described in the continuum limit by the 
Gaussian model with the coupling constant (7 = 1/2 by explicit calculation of the toroidal 
partition function 0. Effect of domain wall interaction has also been studied in the fermion 
model derived from an approximation to the ANNNI model [0 and in a phenomenological 
model 0. 

In this paper, we consider an exactly solvable interacting domain wall model derived 
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from the triangular lattice antiferromagnetic Ising model (TAFIM). It is well known that 
the TAFIM with the only nearest neighbor coupling has infinitely degenerate ground states 
due to frustration on each elementary triangles. Each ground state can be mapped into 
a configuration of covering the plane by three types of diamonds. Blote and Hilhorst ^ 
introduced a solid-on-solid model derived from these configurations. Regarding two types of 
diamonds as domain wall excitations, one also obtains a striped domain wall configuration. 
Blote and Hilhorst utilized this connection to obtain exact solution to the non-interacting 
domain wall problem. As the fugacities of walls change, there is a phase transition from an 
ordered phase to the critically disordered phase which is described by the Gaussian fixed 
point with the coupling constant g = 2. The nature of the transition is found to be that of 
the PT transition . Nienhuis et al. identified various spin wave and vortex operators of 
the Gaussian model in terms of the solid-on-solid model and argued that infinitesimal next 
nearest neighbor (nnn) interactions and magnetic field in TAFIM would change the coupling 
constant g of the Gaussian model. From this they suggested a schematic phase diagram in 
the parameter space composed of the nearest neighbor interactions, the nnn interactions 
and the external magnetic field. More recently the effect of the external magnetic field on g 
has been studied by Blote et al. ||TD[ and the behaviors predicted in Ref. is confirmed. 



We show in Sec. |T|that the ground state configurations of the TAFIM under the general 
boundary conditions are equivalent to the striped domain wall configurations. When the nnn 
interactions in the TAFIM are turned on in an anisotropic manner, they correspond to extra 
energies between adjacent domain walls. Only the same types of walls can touch each other 
and there are two types of wall interactions. We also show in Sec. |I| that the striped domain 
wall configuration is exactly mapped to the arrow configuration of the 5-vertex model. But, 
if both types of wall interactions are present, the Boltzmann weight cannot be represented 
by a product of vertex weights. However, when only one type of domain walls interacts each 
other, it can be written as a product of vertex weights and the partially interacting domain 
wall model reduces to the general 5-vertex model. 



In Sec. |Tl|, we diagonalize the transfer matrix of the 5-vertex model using the Bethe 



Ansatz method. We develop Bethe ansatz solutions both for domain wall and domain wall 
hole. From these solutions, we obtain full phase diagram of the partially interacting model. 
The phase diagram displays the C and IC phases separated by the PT transition and the 
first order transition. It also exhibits a new commensurate phase where the domain wall 
density is locked to the value 1/2 for a range of the chemical potential of the wall. This 
phase does not appear in the non-interacting domain wall models and is a feature resulting 
from the domain wall interactions. This is akin to the antiferromagnetically ordered phase 
of the ANNNI model. 



In Sec. jl^, we investigate the critical properties of the IC phase. It is shown that the 
interaction between domain walls causes a continuous variation of the coupling constant g 
of the Gaussian model resulting in non-universal critical behaviors. It is studied by analytic 
perturbative calculations and numerical calculations. We discuss and summarize our result 
in Sec. and present discussions on the Yang-Baxter equation and the calculation of the 
modular covariant partition functions of the T = TAFIM under the general boundary 
conditions in Appendix ^ and respectively. 



II. TRANSFER MATRIX FORMULATION OF INTERACTING DOMAIN WALL 

MODEL 

We write the hamiltonian 7i including 1/kT of the TAFIM with the nearest and next 
nearest neighbor interaction as 

n = -J2iK + ^a)s^S, - '^aS^Sj - KAf (l) 

<ij> 

where Sj = ±1 is an Ising spin variable at site i, the first (second) sum is over the nearest 
(next nearest) neighbor pairs of sites, K + 6a {sa), a = 1,2,3, are the anisotropic nearest 
(next nearest) neighbor couplings whose index a depends on the direction of the bond <ij > 
») as shown in Fig. 1(a), and finally A/" is the number of lattice sites. 
Monte Carlo simulation and other studies []n],|T2[ show that this system has rich critical 
phenomena in the full parameter space. But, we will only consider the zero temperature limit 
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of this system. By the zero temperature hmit, we actually mean the infinite coupling limit 
K —oo leaving 5j's and e^'s finite. Eq. (|1]) in this limit will be called the T = TAFIM. 
Here, only those configurations which have precisely one pair of parallel spins around each 
elementary triangle are energetically allowed. Though this imposes much restriction on the 
spin configurations, it is important to study this limiting case because the T = TAFIM is 
equivalent to many interesting problems, e.g. diamond and/or dimer covering problem |^ 
and triangular solid-on-solid model |]9|JT3[]. Moreover, the T 7^ behavior of the TAFIM can 
be inferred from the T = behavior. 

Here, we will show that the T = TAFIM with nnn interaction is equivalent to the 
interacting striped domain wall model where the nnn interaction Sj (j = 1, 2) plays the role 
of wall-wall interactions. If we draw lines between all nearest neighbor pairs of antiparallel 
spins for a given ground state configuration of the TAFIM, the resulting configuration is that 
of a covering of the plane by diamonds. Fig. 2 shows a typical TAFIM ground state and 
its corresponding diamond covering configuration. The three types of diamonds are called 
as type 1,2 and 3, respectively as shown in Fig. 1(b). Strictly speaking, there is two-to-one 
correspondence because of the global spin reversal symmetry of the TAFIM in the absence 
of magnetic field. 

From a diamond covering configuration, a striped domain wall configuration is obtained 
by regarding the diamonds of type 1 and 2 as domain wall excitations. Type 3 diamonds 
are regarded as the vacuum. Thick lines on the faces of type 1 and 2 diamonds in Fig 1(b) 
and Fig. 2 visualize the domain walls. A section of domain walls which is obtained from the 
diamond of type 1 and 2 will be called the domain wall of type 1 and 2, respectively. Two 
walls are defined to be interacting when their sides touch each other. The thermodynamic 
parameters which control the equilibrium property of interacting striped domain wall system 
are the fugacities Xi and X2 of domain walls of type 1 and 2, respectively, and the fugacities 
Hi and ?/2 for each pair of adjacent domain walls of type 1 and 2, respectively. Note that 
different types of domain walls cannot be adjacent. The partition function Z^^ for the 
interacting domain wall model is 
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Z,,^=Y^x\-xl-y'^y':i (2) 

where the summation is taken over all striped domain wall configurations and is the total 
length of domain wall of type i and is the total number of incidents where domain walls of 
type i touch each other and share a side, i.e. the number of wall-wall interactions of type i. 

When £3 = in Eq. (P, the energy of the T = TAFIM can be written in terms of nj 
and li. Nearest neighbor interactions contribute simply 

^ {-bi + bj + bk)ni 

j=l,2,3 

where (z, j, fc) is the cyclic permutation of (1, 2, 3) and 723 = M — n\—n2 is the total number 
of type 3 diamonds. From now on we set ^3 = without loss of generality. To relate the 
nnn interaction energies to /j, consider first the bonds connecting nnn pair of sites along the 
direction 1. (See Fig. 1(a).) They cross either (a) two type 2 domain walls or (b) two type 
3 diamonds or (c) one type 2 wall and one type 3 diamond or (d) one type 1 domain wall. 
These possibilities are shown in Fig. 3. If we let na,nh,nc and rid be the number of cases 
(a), (b), (c) and (d), respectively, the bonds contribute ei{na + rib — ric — n^) to the energy. 
But one can easily identify = ni and Ua = h- Moreover each of type 2 wall is crossed by 
two nnn bonds so that it appears twice in the list of Fig. 3 while the total number of type 
2 walls counted in Fig. 3 is 2na + n^. Thus 2n2 = 2na + rtc. These relations, together with 
the sum rule Ua + rih + ric + Ud = J\f, give the energy 

ei(Ar + 4/2 -4^2 - 2ni) . 

Similar counting holds for nnn bonds along the direction 2. 

Putting these together, the ground state energy of Eq. (|l|) for £3 = becomes 

Eo = -2(5i +ei + 2e2)ni - 2(^2 + £2 + 2^1)^2 

+Ae2li + Aeik + AfiSi + 52 + ^1+ £2) • (3) 

Thus the fugacities for the interacting domain wall model are related to anisotropic coupling 
energies of the TAFIM model as 
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xi = exp [2{6i - S3) + 2ei + 452] 
X2 = exp [2(^2 - 63) + 262 + 4£i] 

(4) 

yi = exp [-452] 
1/2 = exp[-4ei] 

Next, we show that to each striped domain wall configuration, one can assign a vertex 
configuration. To do this we deform the triangular lattice into the square one as shown in 
Fig. 4. One then finds that there are 5 types of unit squares. Fig. 5 shows them together 
with assignment of vertex configurations. If one works under the ice rule, the assignment of 
vertices shown in Fig. 5 is unique modulo the arrow reversal. 

In this way, we obtain one-to-one correspondence between striped domain wall configu- 
rations and bond arrow configurations satisfying the ice rule. Vertical up arrow indicates 
the presence of a domain wall. In the TAFIM language, vertical up arrows correspond to 
horizontal nearest neighbor spin pairs which have opposite signs and right arrows correspond 
to vertical nearest neighbor spin pairs which have the same signs. The absence of the third 
vertex (or the fourth upon arrow reversal) is a result of our deforming the triangular lattice 
in the manner shown in Fig. 4. If it were deformed in the opposite direction, it is the first 
vertex (or the second upon the arrow reversal) which does not appear. In any case, one 
obtains the 5-vertex model configurations. 

The 5-vertex model on the square lattice is obtained from the 6-vertex model by sup- 
pressing one of the first four vertices. The 5-vertex model with special choice of its vertex 



weights was first considered by Wu |]T4[ as a limiting case of the 6-vertex model and is 
studied in connection with the non-intersecting directed random walk [|1^ and the directed 
percolation problem in three dimension [|T6|. Recently, Gulacsi et al. [|T^ studied its phase 
diagram for a special case. The general 5-vertex model is obtained by assigning arbitrary 
vertex weights to each type of vertices but there are only three independent parameters 
since the vertices 5 and 6 always occur in pairs along a row under the periodic boundary 
condition and a global rescaling of weights introduces only a trivial factor. 



The partition function Z^_^ of the 5-vertex model is 

where the summation is taken over all arrow configurations and Ni is the number of the z-th 
vertex appearing in an arrow configuration. Unfortunately, if we assign Boltzmann weight 
of an arrow configuration as a product of local vertex weights, we can not treat the fully 
interacting domain wall model [e\ 7^ 0,52 7^ case). However, if we restrict ourselves to the 
special case ei = (7/2 = 1), then Eq. (g) can be expressed in the form of Eq. (^. From now 
on, we consider the partially interacting model where only domain walls of type 1 interact. 

The vertex of type 1 represents the case where the two domain walls of type 1 are 
interacting. So, if we choose vertex weights as 





= XiVx 




W2 


= 1 


(6) 






= X2 




W5 = Wfs 


= 





the partition function Z^_^ becomes the same as that of the partially interacting domain wall 

system; 

z,_^ = J2xTxTy[' . (7) 

We study the 5-vertex model using the transfer matrix. Suppose the lattice has M 
rows and columns, and periodic boundary conditions are imposed in both directions. Let 
at = (ai, ■ ■ ■ ,a7v) denote the state of vertical arrows of one row. Then, as usual, we can 
write the partition function Zg.^ as 

Z,. = Tr Tf, (8) 
where Tg.^ is the 2^ by 2^ transfer matrix with elements 

TV 

Tg..(a,/3)= nW^(/^*'«*IA,/i^+i) ■ (9) 

{,i,=±i}i=i 
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In Eq. (I), 1/) is the weight of the vertex configuration in the standard nota- 

tion 10 . Let Tl (Tr) be the transfer matrix of the 5-vertex model with the first horizontal 
arrow fixed to the left (right). This can be written graphically as 



/3i 132 



13. 



'N 



Tl{<x I (3) 



I3l P2 



13^ 



'N 



Then the transfer matrix can be written as 



■R 



(10) 



From the ice rule, the number of up arrows on a row and right arrows on a column are 
conserved. In the language of domain wall, the number of up arrows per row corresponds to 
the number of domain walls per row and the number of right arrows per column corresponds 
to the number of type 1 domain walls per column. We will call them Q and fi, respectively. 
From the conservation of Q, Tj.^ is a direct sum of submatrices labeled by Q which only act 
on the subspace with Q domain walls. Thus, 



(Tfi,Q 



■L,Q) 



(11) 



where stands for the direct sum and T/j q {Tl,q) denotes the sector Q of T/j (T^). 

The partition function Z^_^ of the 5-vertex model is obtained from the partition function 
^TAFiM of the T = TAFIM as follow. Suppose the triangular lattice has M rows and N 
columns as in Fig. 2 under the boundary condition (yU, u) defined by 



Si,M+l 



(/i = 0,l) 
(z/ = 0,l) 



(12) 



(-1)^5^,1 

where /i, u are (1) for periodic (anti-periodic) boundary condition. Let s = (si, ■ ■ ■ , sn) 
denote the spin state of one row. It can also be represented by (si, a.) where ck = (ai, ■ ■ ■ , otv) 
and ai = — SjSj+i. With the identification of a = 1 (—1) to the up (down) arrow in the 
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i-th vertical bond in the dual lattice, one notes that a is the arrow configuration of a row 
of vertical bonds of the corresponding 5-vertex model configuration. The transfer matrix 
Ttafim of the T = TAFIM is defined through its matrix element Ttafim(si,q: | s\,a') 
which is the Boltzmann weight for two successive row configurations {si, a) and {s'l, a') with 
the boundary condition u along the horizontal direction. Due to the global spin reversal 
symmetry, it takes the block form 

^ + T+_ ^ 
^ T++ - ^ 



-•- ++ -•- +- 



(13) 



where T^^' (s, s' = ±) is the matrix whose elements are Ttafim (-5, ct \ s', a') and ~ denotes 
an equivalence up to the similarity transformation. We use the fact that Tgs' = T^s-s'- If 
we denote the partition function of the T = TAFIM under the boundary condition (/i, u) 
as -Ztafim, it can be written as 

2;M = Tr R'' [t^*"^ 1^ (U) 

^TAFIM [-^TAFImJ V-^^y 

where R is the spin reversal operator. 

Since the sign of spin reverses by crossing each domain wall in the horizontal direction, 
spin configurations under the boundary condition u = (1) yield only domain wall configu- 
rations with Q even (odd). Therefore, T^.^ and T+_ in Eq. (p!3D are Tji and T^, respectively, 
of the 5-vertex model restricted to Q even (odd) sectors for u = (1). This shows that the 
transfer matrix spectra of the two models are not identical. Only the even or odd Q sector 
of Eq. (|TT]) are present in Eq. (|I^ while the latter includes the block T_|_+ — T+_ which is not 
present in the 5-vertex model. And the spin configurations under the boundary condition 
/i = (1) yield only domain wall configurations with Q even (odd) assuming that M is even, 
since the sign of spin changes in every step except when crossing a type 1 domain wall in 
the vertical direction. If M is odd, even (odd) fl corresponds to = 1 (0). So, the partition 
function Zr^_^ of the 5-vertex model is given by 

7 =if?'(0'0) -L T-CO,!) , y(0,l) , ;r(l,l) A c^ 

2 V ''^^^'^ TAFIM ~ ^TAFIM ~ '^TAFIM J 
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where the factor 1/2 accounts for the two-to-one correspondence. This relation will be used 
in Sec. |V| to obtain the toroidal partition function Z^_^ of the 5-vertex model. 

III. PHASE DIAGRAM 

The five-vertex model transfer matrix can be diagonalized by the Bethe ansatz method 
as a special case of the general six-vertex model [|I^. Its phase diagram has recently been 
calculated by Gulacsi et al. |[T^ for the special case of wi = W2 Q- In this section, we 



generalize it to the full three dimensional parameter space and also calculate the domain 
wall densities. We also discuss types of solutions of the Bethe ansatz equation (BAE) of the 
5-vertex model. 

The eigenvalues of the transfer matrix Eq. (|]) in the sector Q (7^ A^) are given by 0Jlj 



aq = wr^w^ n 1 + (16) 



where the set {zi, Z2, ■ ■ ■ , zq} are the solutions of the BAE 

1 - Azi 



= (-i)^-^nY^ ' ^■ = i'2'---'^ (17) 



with 



W2W4 

All z/s should be distinct. When Q = N, An = Wi + . 

An alternative expression for the eigenvalue which is useful for Q > N/2 is given by 

^Q = ^illi \+w^[[\W2 —\ (19) 

where Q = N — Q is the number of domain wall holes and the set {zi, Z2, . . . , zq} is again 
given by Eq. (|1^ with Q replaced by Q. We call Eq. (|1^) and Eq. (^) the domain wall 



^Gulacsi et al. use the notation W2 = 0- Their work and ours are related by the transformation 
wi W4, W2 ^ ws and W5 wq. 
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representation and the domain wall hole representation, respectively. Using Eq. into 
Eqs. ([1§), QTSP and (|19|) gives Ag in terms of the domain wall parameters as 

Q 

AQ = x^l[{l + azj) (16') 
j=i 

or 

A«^.?(A..)4{^^^)..«ft{l-^^^) (1.) 
where a is the ratio of two domain wall fugacities 

a = xi/x2 (20) 

and 

A = xi(yi - 1)/X2 = a(i/i - 1) . (18') 



We will call A the interaction parameter. It is positive for attractive interaction and negative 
for repulsive interaction between domain walls. 

Defining the momenta {pj} by zj = e*^^, Eq. ( p!7| ) also takes the familiar form 

Q 

Npj = 27rl,+Y.e{p,,pi) (21) 
1=1 

where 

1 _ Ap*P 

gie{p,,) ^ _ ^ (22) 
1 - Ae*? ^ ^ 

and Jj's are half-integers for even Q and integers for odd Q ranging from — (A^ — l)/2 to 
(A^ — l)/2. Different eigenvalues come from different choices of the set {Ij}- 
The BAE may take another form. If we define 

^ = ^EMi-A^O , (23) 

then the BAE becomes 
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= (-1) (1 - Az,y e-'^' (24) 

where q = Q/N is the domain wall density. This equation gives zj^s as a function of s which 
should, in turn, be determined from its defining equation 

Note that the BAE (Eq. (|l^) arises from the periodic boundary condition on the wave 
function of T^_^ . It is also interesting to consider another boundary condition, say, the 
anti-periodic boundary condition. The effect of the boundary condition is to shift domain 
walls out of the A^-th site to the first site with appropriate phase factor 1 (—1) for periodic 
(anti-periodic) boundary condition. The shift operation is done by the operator T/j. So, if 
we impose anti-periodic boundary condition, the resulting matrix we diagonalize is T^ — T^^. 
In this case, the expression for eigenvalues remains the same except the fact that /-,• should 
be integers for even Q and half-integers for odd Q. So, we can obtain whole spectrum of 
the transfer matrix of the T = TAFIM from the transfer matrix of the 5-vertex model 
under periodic and anti-periodic boundary conditions. Note that the anti-periodic boundary 
condition here is different from that which reverses the sense of horizontal arrows. 

The free energy in the language of the domain wall physics is a function of xi, X2 and yi 
through Eq. (|^). From now on, we regard it as a function of A, X2 and a = Xxjxi- Since the 
free energy is given by the maximum eigenvalue of the transfer matrix, /(x2, a. A), the free 
energy per site in the units of fcT, is written in the form 



/(a;2,a. A) = - lim max 



.A^ y 



= — max [glnx2 + K(g)] (25) 

where ^(g), which will be called the configurational free energy, is given by 

1 



Kiq) = lim max 

N-*oo {zj} 



(26) 



Here, {zj}'s are the solutions of the BAE. The equation of state which relates the equilibrium 
domain wall density g as a function of thermodynamic parameters is given by the relation 

q(x2,a,A) = ^ (27) 
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where Qo is the value of Q at which Ag attains the maximum value. The equation of state 
can be rewritten as 

d 

— K{q, a, A) = - In X2 (28) 

if K,{q) is a differentiable and convex function. The configurational free energy K,{q) is a 
Legendre transformation of /. That is, it is a free energy as a function of domain wall 
density while / is a free energy as a function of the domain wall fugacity. 

We now classify types of solutions of the BAE corresponding to the maximum eigenvalue. 
First, consider the case — 1 < A < 1. This region contains the non- interacting case with 



A = which is considered in |14|. In this case, the G function defined in Eq. (|22|) is identically 
0. So, any set {Ij} of Q different numbers are solutions of the BAE and the solution giving 
the maximum value of Ag is {Ij} = {-{Q - l)/2, -{Q - l)/2 + 1, . . . , (Q - l)/2}. We 
assume that this set {Ij} still gives the maximum eigenvalue even after turning-on of weak 
interaction and remains so in the whole region — 1 < A < 1. This assumption is tested by 
direct numerical diagonalization of the transfer matrix with N up to 15. We call this type 
of solution as the free magnon type. 

When |A| > 1, there appear other types of solutions. Assume that the solution is of the 
form, 

Zj = A'^Zj , j = l,...,iV+ 



^3 ^3 



j = (Ar++l),...,(Ar++iVo) (29) 



z. 



h = X-TTT^ , j = {N^+N, + l),...,Q 



A Ai+'*^- ' 

where aj and bj are constants greater than and zj^s are assumed to remain of order 1 
as |A| oo. In other words, of Q Zj's, are diverging, A^^o remain finite and A^_ = 
Q—{N^+Nq) vanish inversely as |A| —>■ oo. Then, the necessary condition that this set 
should be a solution of the BAE is either (i) A^+ = A^_ = so that all z/s are of order 1 
or (ii) No = 0, bj = {Q - N+)/N+ and = {N - Q)/N+. We also call the first type as the 
free magnon type while the second type will be called as the hounded magnon type. Of these 
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possibilities, one can easily show that the configurational free energy is realized by the free 
magnon type if A ^ — oo and the bounded magnon type with A^_|_ = 1 if A ^ oo. We find 
numerically this feature also persists for all A in the range |A| > 1. 

It is very difficult to obtain the equation of state (Eq. (pH])) analytically for whole range 
of parameters X2,a and A. But, the phase boundaries which separate the commensurate 
phases with domain wall density or 1 from the incommensurate phase can be obtained 
if we solve the BAE in the g — > or 1 limit. Apart from the C phases with domain wall 
density or 1, there appears new a C phase with q = 1/2 if a is large so that A can take 
values less than some critical value Ac- (Below, we will see that Ac takes the value —4.) 
Consider the case where domain walls of type 1 are much more favorable to form than those 
of type 2 (a = Xi/x2 ^ 1) and there are repulsive interaction between them (jji < 1) so that 
the interaction parameter A is less than A^. Then, the most energetically favorable state 
for g < 1/2 is the state where there are only type 1 domain walls with no adjacent pairs. 
But, if q is larger than 1/2, there should appear type 2 walls and adjacent pairs of domain 
walls of type 1 whose energy costs are large. So, it is expected that there is a discontinuity 
in X2 which controls the total number of domain walls across the g = 1/2 line. The phase 
boundary of the g = 1/2 C phase can be also determined analytically. Our results for the 
phase boundary are given by Eqs. (^), (il), (|38|) , (^6|) , ( ^6|) and (|67|) and are illustrated in 



Fig. 7. 

(i) A < 1, g = 

In the region A < 1, the free magnon type solution in the g — limit is 

\- 0{q') (30) 



where 6j = 27i:Ij/N and {Ij} is a set of integers or half-integers depending on the parity of 
the domain wall number Q. The maximum value of K,{q) is obtained if we choose the set 
{Ij} = {-{Q - l)/2, -(Q - l)/2 + 1, . . . , (Q - l)/2} and the next largest values of fi:(g) are 
obtained by using the set {/j = Ij + m} which is a shift of the set {Ij} by an integer m. 
With this solution, the configurational free energy k is given from Eq. (|26D by 
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= q ln(l + a) 



/ In 


1 + ae'^ ( 


1 — nq 





avr 



1 — X 



(31) 



6(1 + a) 

This, together with Eq. (p8|), imphes that g = if X2 < 1/(1 + a). Thus, we obtain the 
phase boundary X2 = xqc between the q = C phase and the IC phase as 

1 



1 + a 



(32) 



or equivalently Xi + X2 = I. For X2 shghtly larger that Xqc, Eq. (|3l|) gives 



'2(1 + a) 



(lnx2 - Inxoc) 



1/2 



(33) 



The domain wall density thus shows the square root dependence on domain wall formation 
energy which is the general character of the PT transition. This type of singularity is 
originated from the fact that the leading contribution K,{q) aside from the linear term is of 
order q^. It is originated from the entropy reduction due to the collision of domain walls pO . 
(ii) A < 1, g = 1 

Next, consider the case near q = 1. In this case, it is easier to consider the BAE for 
domain wall hole rather than domain wall. Inserting Eq. (|30| ) with q replaced by q into 
Eq. ( [19'D , we obtain K{q) near q = 1. g = l — gisa domain wall hole density. There are two 
cases to consider depending on whether A + a> lorA + a< 1. When A + a > 1 , the first 
term in the right hand side of Eq. (|19'|) dominates and hence 

Kiq) = q ln( A + a) + 1 g In [ (^^^^ 



, ^ , (A + a)(A + a-l) f n'^(A + a) ,^ 
de = ln(A + a) - g In ^ '-^ - ——-^ ^ + 0{q^) 



(34) 



a 6 (A + a-l)2 

From this, the phase boundary X2 = xic between g = 1 C phase and IC phase is given by 



(A + a > 1) 



(A + a)(A + a- 1) 
and the equilibrium domain wall density near xic is given by 



(35) 
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g ~ 1 — 



2(A + a-l) 



(Inxic - 111x2 



,1/2 



\ 7r2(A + a) 

Similarly, when A + a < 1, the configurational free energy is given by 



(36) 



Q 



K 



1 - 



a 



(A + a)^,--l 



1-A 1 3 a7r2(l - aA - A2) 



(37) 



1-A- a 6^ (A + a- 1)2(A- 1)2 

from which the phase boundary X2 — Xic and the equilibrium domain wall density near xic 
are given by 

1-A 



xic = 



1- A-a 



(A + a < 1) 



(38) 



and 



g ~ 1 — 



^^^:r^r^\'!^ in..-in..)V2 



(39) 



\\ an^il -aA- A^) 

So, we conclude that when A < 1, there are commensurate phases with domain wall density 
for X2 < xqc and domain wall density 1 for X2 > xic- In between, the equilibrium domain 
wall density increase smoothly as X2 increases as long as A > —4. Fig. 6 shows a typical X2 
dependence of q for the case of a = 1. The curves are obtained numerically by solving the 
BAE for N up to 150. The case A < — 4 will be considered later. 
(iii) A > 1 

Now, consider the case where there is a strong attractive interaction between domain 
walls so that A > 1. The solution of the BAE maximizing k is of the bounded magnon type 
with — 1. The exact solution of the BAE is easily obtained from the transformation of 
z to z defined by 

zi = ^iA(«-i) 



Then the BAE for Zj becomes 
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For macroscopic number of and Q, the values of Zj/A'^ ^ and A are exponentially 
small and may be neglected. Thus the solution is 



(41) 

where 9 can take the value from the set ^ x {1, . . . , A^}. The corresponding right-eigenvector 
I A) of the Tr^q + Tl,q is 

\\) = C e"^''A-S^=2("'-"iVi>--->«Q) (42) 

ni<--<nQ 

where C is a normalizing constant. Here, {nj}'s denote the position of up arrows or equiv- 
alently domain walls. It is obvious that these states represent bounded domain wall states 
because the components of eigenvector decay exponentially as the distance between domain 
walls becomes large. In fact, one can calculate the mean distance of the last domain wall 
from the first one. The mean distance {uq — ni) of Q domain wall system is given by 

<"-'-"->= ^"%r"'^ (43) 

where fij is the position operator of the j-th domain wall. Inserting the eigenket to the above 
expression and after some algebra, we find that {nq — rii) is equal to Q for macroscopic 
number of domain walls. Thus, we can interpret this state as the bounded domain wall 
state. 

This solution yields the exact configurational free energy which is obtained from the 
solution Eq. (|l]) with ^ = 0; 

K{q) = q\n{A + a) . (44) 
And the free energy / for A > 1 is simply 



— Pf = ^^^^ [l In X2 + g ln( A + a)] . (45) 

The maximum value is obtained at g = if 0:2 is less than 1/ ( A + a) and at g = 1 if 
X2 is greater than 1/ ( A + a) . So there is a first order phase transition between the two 
commensurate phases when X2 is at the critical fugacity Xc, where 

Note that the condition wi = W2 used in is amount to the condition X2 = 1/(A + a) so 



that the first order transition for A > 1 could not to be seen in |T^. We have thus found 
the phase boundary of the C phase with domain wall density q = and 1 and the nature 
of the phase transition. We present the resulting phase diagram in Fig. 7(a) for the case of 
a = 1. 

(iv) A < -4, g = 1/2- 

As discussed before, we expect that / has a singularity in X2 at g = 1/2 if A is large and 
negative. To see q dependence of X2 near q = 1/2, we should evaluate the configurational 
free energy K,{q) near g = 1/2. Gulacsi et al. |]l^ used the root density function p{p) to 
find the g = 1/2 phase boundary when a;2(A + a) = 1. We employ the same method to the 
general case. 

p{p) is defined so as Np{p)dp to be the number of the roots of the BAE (Eq. (|17D) with 
z = e*^ in the interval {p, p + dp) in the complex p plane. We stress here that the roots do 
not lie on a straight line in the complex p plane. For domain wall density g, p{p) is given 

by m 

In Fig. 8, we give a typical root distribution of the BAE in the complex a plane which is 
related to p as 

= 1 - Ae'P . (48) 

The root density function p{a) in the a plane is given by 
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and the variables A, B and D in Fig. 8 are given by 

/(vr + A) 
[infe^-^-^ 



g = Im 

1 - e'^-^ 



Re 



In- 



A 

D = 1 - AC 



+ -Re 

TT 



1) 



+ 



AB 



TT TT 



— — — sin = 



(50) 



n=l 



where C is determined from the equation 

(1- AC)" 



C 



exp 



^E7=i ln(l- Ae^f. 



(51) 



Near g = 1/2, they take the values 

/I = TT + (g - m^\^ + 0{{q - 1/2)2) 

B = lnA + 0((g- 1/2)2) 
lA ' 



(52) 



D = 1 
with 



1 



'1 



4 
AJ 



(g-l/2)ln(dA)/ 



d-1 2d 



+ 0((g- 1/2)2) 



A^-(|A| 



(|A| -2)2-4 



(53) 



With this knowledge, we can calculate K,{q) near q = 1/2. First consider the case g < 1/2 
where ^(g) is evaluated from Eq. ( [L6'|) . 

Q 

N 

'a + A — ae*" 



ae 



da p{a) In 



(54) 



ic ■ ■ ■ \ A ; 

where the integration should be taken along the contour C shown in Fig. 8. But, the contour 
can be deformed to the straight line C since the integrand is analytic in the shaded region. 
Then ^(g), up to the first order in (g — 1/2), is given by 



K{q) 



1 

271 J -A 



-q + 



e^VA - 1 



In 



'a + A-ae'^/y 



-lna+(g-l/2) 



a , / a + A , 
In + 2 In 1 + A 



lAIA 



a 



0{{q -l/2f 



(55) 
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Note that /t(l/2) = |lna independent of A. This imphes that there are only type 1 domain 
walls without adjacent pairs of them at q = 1/2 phase. From the Eq. (^), we see that 
g = 1/2 phase starts at X2 = a;_ where 



In a;_ = — In -— - — 2 In 
AA 



1 + ^A- 



(56) 



(v) A < -4, g = l/2^ 



Next consider the case q > 1/2, where the configurational free energy is evaluated from 
the Eq. (|19]). For convenience, we define two quantities A and B; 

Q 



B^j^j:in{l-{a + A)e^^^). (57) 

Then the configurational free energy is written as 

K{q) = max{Ki(g), /tij(g)} (58) 

where = A — B and = lna + gln((a + A)/a) — ^B. The quantity A and B can be written 
as a contour integration in the complex a plane in the thermodynamic limit N —>■ oo. 



A = J da ia p{a) 



B = ldapia) ln(^^^^^^). (59) 

Since the integrand in A is analytic in the shaded region, the contour can be deformed to 
the straight line C and the integration results in 

^=iln|A|-(g-i)lnA + 0((g- 1/2)2) _ (gQ) 

In other to calculate B, there are three possible cases to consider; 
(a) 1/D < A < (a + A)/a 

The integrand is also analytic in the shaded region and the contour C is replaced by the 
straight line C It yields 
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Ba = ^ Ha + A) + (g - 1/2) 



, a + A , / 
In + 2 In 1 



aX 



O Uq- 1/2)'' 



(61) 



|A|A V a + A^ 

(b) l/D < (a + A) /a < A 

In this case, the branch cut intrudes into the shaded area. Therefore, upon changing 
C to the straight hne, one need to subtract the contribution around the branch cut. The 
results is 



Bb = ^\n{a + A) + {q- 1/2) 



'a + A / aX 

In + 2 In 1 + 



+ 0{{q- 1/2)'' 



(62) 



|A|A \ a + A^ 

(c) (a + A)/a < l/D <X 

In this case, the contour can be deformed to the straight line C as in the case (a) and 
the integration results in 

a + A' 



i3c = ^ln^ + (g-l/2) 



In — + 2 In 1 + ^ 
A V aX 



+ (g-1/2)' 



(63) 



For each case, the quantity and kl take the following values at g = 1/2; 



(a) = l\n^ , Kn 



I In a 



(b) 



In- 



2 a+A ' '^R 2 

(c) = ilna , 



i Ina 



(64) 



llna + iln^ 



When a + A > 1, only the cases (a) and (b) in Eq. ( |64D occur and > hl sX q = 1/2. 
So near q = 1/2, the configurational free energy K{q) is determined by KR{q) and is given by 



«:(g) = ^lna+(g-l/2) 



, Xa , / a + A 

In — + 2 In 1 + — 

A V aA 



+ 0{iq-l/2)' 



(65) 



When a + A < 1, the cases (b) and (c) occur. One can easily find that if (a + A)/a < l/D 
then kl > hr and if (a + A)/a > l/D then kl < i^r- So, near g = 1/2 the configurational 
free energy is 

KL{q) , if (a + A)/a < l/D 



K{q) 



KR{q) , if (a + A)/a > l/D 



-lna + (g-l/2) 



, Aa , / a + A 

In — + 2 In 1 + — 

A V aA 



0{{q-l/2f 



(66) 
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From the Eq. ( |65D and (|66D , we see that the q = 1/2 phase ends at X2 = x+ where 

lna;+ = -ln^-21n(l + ^) (67) 

X- (Eq. (|56D) and a;+ (Eq. (|67D) defines the phase boundary between the g = 1/2 C 
phase and the IC phase and the domain wall density is locked at q = 1/2 for the range 

x_ < X2 < x+ . (68) 

Since x_ and x+ merge at = —4, this phase appears only when a > 4. Fig. 7(b) shows 
the full phase boundaries for a = 7 and Fig. 9 shows the domain wall density as a function 
of lna;2 for A = —4 and —5. 



IV. THE CRITICAL PROPERTIES OF THE IC PHASE 

The conformal field theory predicts that the operator content of a critical phase is related 
to the finite size correction to the eigenvalue spectra of the transfer matrix When we 
write an eigenvalue Aq, of the transfer matrix for a lattice of width N as e~^° , then takes 
the form at the criticality, 

= Ar/^ + ^ (a„ + A, -^)Csin0 + ^(A, - A,)Ccose + o(^^) , (69) 

where c is the central charge, (^Aq,, A^^ are the conformal dimensions of the operator corre- 
sponding to the a-th energy eigenstate, ( is the anisotropy factor, 9 is the anisotropy angle 
and finally foo is the non-universal bulk free energy per site in units of kT [pl| , p^ . 

The toroidal partition function (TPF) Z is defined as the order 1 part of the partition 
function Z for conformally invariant system of N columns and M rows. It follows from 
Eq. (§yp that 

Z= hm ye-^^^Ve"'^"''^- 

M/N=fixed ° 

V- r 27rMC 

= Eexp — ^ 

a L 



|(^A^ + A„ - sin^ + 2 (A„ - A^) cos^ 
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where g, the modular parameter, is given by 

q = e^^'^ (71) 

with 

r = |Ce^(-^) , (72) 

q is the complex conjugate of q and the sum is over the infinite set of levels whose energy 
Ea scales as Eq. (|69|) . In the first part of this section we use the notation q to denote the 
modular parameter (Eq. ([7TD). This is not to be confused with the domain wall density. 

For the Gaussian model compactified on a circle, or equivalently, the symmetric six- 
vertex model in the continuum limit, the TPF under periodic boundary conditions in both 
directions is given by the c = 1 Coulombic partition function 

^ci9) = Y^ E (73) 

1 / n Y 

where g is the so called Gaussian coupling constant, A^^ = h ^/gm and 77(g) is 

4 \^ J 

the Dedekind eta function; 

00 

^(g)=g2^n(l-g") . (74) 
n=l 

One can impose U{1) boundary conditions on the six-vertex model instead of periodic 
boundary conditions. In the Pauli spin representation, the twisted boundary condition is 

(^^+1 ± ^^'N+l) = e^^^ K ± ^^^) (75) 
where ip is the twisting angle. The Coulombic toroidal partition function is then modified 

to [g 

Zc{g) = - E e~*"'^'g^"''"-'^/2,r(9)gA„,_(„_^/2^){5) ^^g^ 

\v{'l)\'^ n,rn£Z 

where cp and ip' are the twisting angles in the space and time directions, respectively. 
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After this short review, we now turn to the critical properties of the IC phase. It is 
generally known that the striped IC phase is critical and described by the c = 1 confor- 
mal field theory in the continuum limit In the fermion model approach, Park and 

Widom 1^ calculated exact toroidal partition function explicitly for the free fermion, i.e. 
non-interacting domain wall system and showed that it is of the form of Eq. ( |7BD where 
g = 1/2, (fi' = and (f/2TT is the number of the domain walls per row (mod 1). Note that 
the twisted boundary condition used in Ref. P] has no direct physical meaning. 

For the T = TAFIM without the second neighbor interaction, the central charge 
and the scaling dimensions of several operators are calculated analytically pSf. Since all 



the transfer matrix spectra are known from the Onsager solution in this case, one may go 
one step further and calculate the toroidal partition function explicitly. We present the 
calculation in Appendix |B[ 

When A = 0, the TPF Zr^_^ of the 5- vertex model can be obtained from the 2^tafim of 
Appendix O, by using the relation Eq. (|T5|). The result is 



= ^_ V e-2^-*^™g(-+^)V2g~(-^)V2 ^ 0) . (77) 

This takes the final form after the modular transformation f — > r = — 1/f; 

Z,_^ = e"*"^i'"g^™'"-«o(^=^/^)g^™'-("-«o)(3=i/2) = 0) (78) 

where Qq and Qi are given in Eq. (P28|) and (P29|) , respectively. This is the exactly Coulom- 
bic partition function with the twisting angle (f = 2ttQq and (f' = ttQi. 

Note that this can be also obtained by replacing (m,n) in Eq. ( P30|) by (2m,n/2). In 
Sec. |T|, we gave the relation between the transfer matrices of the T = TAFIM and the 
5-vertex model. Tg_^ contains odd Q sectors while Ttafim contains the spin reversal odd 
sector, Tl — T^. So, one expects Z^,^ can be obtained from ^tafim by adding terms coming 
from the odd Q sectors and eliminating the terms originated from — T/j. Our results 
show that this is exactly done by a simple substitution of (m, n) by {2m,n/2). 
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Eq. ([78|) implies that the IC phase of the non-interacting domain wall model is in the 
universality class of the Gaussian model with coupling constant g = 1/2 regardless of the 
anisotropies in the fugacity of the domain walls. This result is in accord with previous works 
but it confirms the universality in the strongest sense. 

We assume that the effect of domain wall interactions preserves the c = 1 nature through- 
out the IC phase even though it may change the modular parameters, the coupling constant 
etc. Since the coupling constant g determines the critical exponents, its possible dependence 
on interactions over the IC phase is of interest. If we denote the eigenvalue of Tg.^ corre- 
sponding to the m-th spin wave operator in the sector Q by e"'^'"-'^, it is expected to take 
the form in the IC phase 

where Qo = is the average number of domain walls per row. Here and below, q denotes 
the domain wall density. We now calculate g perturbatively in the small A limit and 
numerically for a wide range of A. During the perturbative calculation with |A| < 1, we 
will only consider the isotropic case (a = 1) for simplicity. In this case, the eigenvalue e"'^™''^ 
of the transfer matrix with A = is 

Ke{^.„,«(A ^ 0), ^ . «M! _ ^) , (80) 

where c = 1, = |tan(7rg/2) and (7° = | as given in Appendix ^ and the superscripts in 
g^ and C,^ denote the value for non-interacting case. If we insert pj = rij + uj into the BAE 
where nj = Ij + m is the solution of the A = BAE for the m-th excited state in a given 
Q sector, the resulting equation for uj is, up to the first order in A, 



Uj = iq 

where s is 



s + Ae-'i'e' 



(81) 



-A^e--/- 1 + 4- (82) 
■nq \ QN^ ^ ' 
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that is determined from the condition J2j Uj = 0. With this solution {uj}, we can calculate 
the energy shift SEm,Q = -Em,Q(A 7^ 0) — E^^q^/S. = 0) due to the interaction; 



N 



1 _|_ pi{nj+Uj) 

yin ^ . 

1 + e™j- 



E 



_5g«"-j _|_ 

1 + e™j 



A ( — h — sin vro I + A ( 2 sin^ — — vro sin vro ) ( — 
V 2 27r 7 V 2 ^ V VA^ 



(83) 



Using the value of Re{Em,Q} / N at A = 0, we can write down the energy Re{£'m,Q(A)} up 
to the first order in A. 



1 

N 



Re{E„,Q(A)} = lRe{E„,,Q(A = 0)} 



N 



Q-Qo 

N 



1 



Re{6E^^Q} 
A 

TT 



— ( vrg sm vrg + 4 sm — j 



+2tt 



C° A / . 27rg 

— - H 2 sm 7r(7 sm no 

2g0 271 \ 2 ^ ^. 



TT 



/^ + |^(|(Q-Qor + ^-^-c/12^ 



(84) 



The new anisotropy factor C, the Gaussian coupling constant g and the central charge c are 
obtained by the comparing the last two expressions; 



9 = sinvrg)+0(A^ 



(85) 



c 



2 tan 



(1 - 2Ag 



cos 



2 Trq 



+ 0(A2) . 



The result from the first order perturbation calculation shows that the interaction between 
domain walls causes a continuous variation of the coupling constant g so the scaling dimen- 
sions vary continuously as a function of the interaction parameter A. 

For larger values of A, g can be evaluated numerically by the finite size corrections of 
the eigenvalue of the transfer matrix (Eq. (|69|) ). Suppose the model parameters are tuned 
in such a way that Qq = Nq is an integer. That is, we are considering the case of q being 
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integer multiple of From Eq. (0), g and Csin^ can be evaluated if we calculate four 

eigenvalues Em,Q with (m, Q) = (0,Qo), (0,Qo±l) and (l,Qo)- 



N 



TT 



^sin^^ = _^Re{Eo^Q„+i + - 2Eo,Qo}Re{Ei,Q„ - i5o,Qo}/2 ■ (86) 



Necessary Em,Q's are obtained by solving the BAE for up to 150. The coupling constant 
g obtained in this way is shown in Fig. 10 as a function of q for several values of A for a 
particular value of a = |A| + 0.1. Note that the value of g starts from around 1/2 at g = 
and ends at 1/2 at g = 1 and varies smoothly when A > —4. The values g = 1/2 at q = 
is easily understood since the interaction effect will vanish in these limit. So is the case for 
q = 1 and A + a < 1. When A < —4, the value of g approaches 2 as g — 1/2. The fact 
that g = 2 exactly in the g — 1/2 limit can be derived analytically following the procedure 



similar to that used by Gwa and Spohn |26 



The domain wall model is obtained from the TAFIM by neglecting the spin configurations 
in which three spins on any elementary triangle are in the same state. This excitation driven 
by the thermal fluctuation creates or annihilates two domain walls at a time and causes a 
domain wall density dislocation. (See Fig. 11.) When two dislocations of up-triangle and 
down-triangle occur in pair, the density dislocation remains as a local defect. These pair 
excitations are analogous to the vortex and anti-vortex pair excitations in the XY model. The 
scahng dimension for the density dislocation p is Xo,2 = Ao,2 + Ao,_2 since such excitation 
creates or annihilates 2 domain walls. Since Xo,2 = '^g "we see that when g < I (xo.2 < 2) 
the density dislocation is relevant and destroys the criticality of the IC phase. Therefore, if 
dislocations are allowed with finite cost of energy, the IC phase cannot remain critical and 
becomes the disordered fluid phase. On the other hand, when g > 1 (a;o,2 > 2) the density 
dislocation is irrelevant and the criticality of the IC phase survives. At the boundary g = 1, 
the KT transition would occur. Since the non-interacting domain wall system has g = 1/2, 
the critical IC phase cannot survive from the density dislocation. However as seen in Fig. 10, 
g crosses the critical value g = 1 in the region of repulsive (A < 0) interactions. The Dotted 
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line in Fig. 7(b) inside the IC phase denotes the position where g takes the value 1. 

So, we conclude that the IC phase near the g = 1/2 C phase is stable under the density 
dislocation. This shows that there are three phases encountered if we consider the dislocation 
effect. They are long-range ordered g = 1/2 C phase, quasi-long range ordered IC phase and 
the disordered phase. They are separated by the PT transition and the KT transition. It 
also explains the phase diagram of the TAFIM with the isotropic nnn interaction obtained 



by Monte Carlo simulation [11 



V. SUMMARY AND DISCUSSION 

In this work, we have introduced a solvable interacting domain wall model derived from 
the T = TAFIM with anisotropic nearest neighbor and nnn interactions. The model is 
shown to be equivalent to the 5-vertex model and exact phase diagram is obtained in the 
three dimensional parameter space. It shows C phases where the domain wall density is 0, 
1/2 or 1 and the IC phase in between. 

The IC phase is a critical phase described by the Gaussian fixed point. The Gaussian 
coupling constant g which determines the scaling dimensions of operators is a function of 
the model parameters and changes smoothly from 1/2 at g = and g = 1 phase boundaries 
to 2 at the g = 1/2 phase boundary. As the interaction is turned on, it decreases (increase) 
for the attractive (repulsive) interaction. For strong repulsive interactions, there is a region 
with g > 1 m which dislocation is irrelevant. Therefore the scenario proposed by Nienhuis 
et al. for the effect of the isotropic nnn interaction in the T = TAFIM is partly born out 
in this model. 

We also have shown by the explicit calculation of the TPF of the non-interacting T=0 
TAFIM that it renormalizes to the Gaussian fixed point with the coupling constant g=2. 
This is in accord with previous works. But, the transfer matrix spectra of the 5-vertex model 
with A = and that of the non-interacting T=0 TAFIM are different in that some sectors 
present in one are absent in the other. This re-distribution of sectors or operator content 
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changes g of the 5-vertex model to 1/2 when A = 0. The TPF of the 5-vertex model is found 
to take the form of the symmetric six-vertex model with the twisted boundary conditions. 
Fractional part of the number of domain walls across a row and a column determines the 
twisting angles of C/(l) boundary conditions along the space and time directions, respectively. 

The model considered in this work is rather special in that only one type of domain walls 
interact. For the fully interacting case, say yi = 1/2 = y, one need to rely on less accurate 
numerical methods. The effects of the interaction between domain walls in both direction 
on the phase diagram and the critical properties of the IC phase are of interest and left for 
further works. 
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APPENDIX A: 

In this appendix, we discuss the Yang-Baxter equation (YBE) of the 5~-vertex model 
and alternative parametrization from which the corresponding quantum chain hamiltonian 
is derived. 

The YBE for the 5-vertex model is given by 

(1 ^ R) (R' 1) (1 R") = (R" ® 1) (1 ® R') (R ® 1) (Al) 
where 1 is the 2x2 unit matrix, ® denotes the direct product, R is the 4x4 matrix given 

by 
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R 



^wi ^ 

W5 W3 

W4 Wq 

W2 



V 



(A2) 



and finally R' (R") is the same as R with Wi replaced by {w'- ). When w^, 



YBE has a solution provided 



A 



W1W2 — w^We w[w2 — uj[^Wq 



W2W4 



The solution under the normalization W2 



Wo 



W2 = 1 is 



0, the 



(A3) 



w'( = w[/wi 
w'^ = 

w'l = {wiw'^ - W4w[) / (w^we) (A4) 
'^5 = 
K = 

As a result, the transfer matrix of Tg_^ having three independent parameters forms a 
two-parameter family of commuting matrices. Vertex weights of the 5-vertex model used in 
this work is given by Eq. (6). If one parametrizes them alternatively as 



Wi 


= e" 




W2 


= 1 




W3 


= 






= (e^ 


- e")/A 


W5 


= ^6 





(A5) 



and similarly for Wj"s, the transfer matrices T^_^{u,v; A) with different u and v commute, 
i.e. 



[T,_^{u,v;A) , T,_^{u',v';A) ] = 
for all u, u',v, and v'. Eq. ([A4|) with above parametrization becomes 



(A6) 
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w'l = e"-" 

^' = e''(e^'-''-e"'-^)/A (A7) 
w'^ = wl = e("'-")/2 

Standard parametrizations of solutions to the YBE involve the so-called spectral parame- 
ter u with which the YBE displays the difference property; i.e., if R = R(m) and R' = R(m') 
then R" = R(m'— m). At criticality, it gives the physical meaning of the anisotropy angle ||22|| . 
Also, corresponding quantum chain hamiltonian commuting with the transfer matrix is ob- 
tained by the logarithmic derivative at m = 0. We find from Eq. (|A7|) that the 5-vertex 
model also displays the difference property if we set f = 0. This is the special case Wi = W2 
considered in Ref . |0 . 



We calculated the quantum hamiltonian 7i of the one- dimensional quantum spin chain 
by taking the logarithmic derivative of the transfer matrix at m = for the case of f = 0. 
The result is 



n = T7 



du 



u,v=0 

^ r A 1 

= j:{^t s^i + jstsl,\ (A8) 

where Sj is the quantum spin density operator at site i. This non-hermitian hamiltonian is 
similar to the hamiltonian of the XXZ quantum spin chain. The difference is that there is 
no term s~sf_^i in this model. So, there is a net flux of the spin flow from the right to the left 
of the chain. It comes from the anisotropic choice of the vertex weights at the beginning. 

APPENDIX B: 

In this appendix, we present the phase diagram of the T = TAFIM with anisotropic 
nearest neighbor interaction and the toroidal partition functions under the general boundary 
conditions. 

Through the star-triangle relation, the Ising model on the triangular lattice can be 



mapped into the Ising model on the honeycomb lattice ||T8|. Let Ki = K + Si and 
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Li {i = 1,2,3) be the interaction strength (including the factor —1/kT) of the Ising model 
on the triangular and honeycomb lattice, respectively. Then, the partition functions ^tafim 
on a triangular lattice with A/" sites and on a honeycomb lattice with 2A/' sites are related 
by 



(Bl) 



provided K and L satisfy the star-triangle relation: 

exp [K1S2S3 + K2S3S1 + K3S1S2] = -R ^ exp [t{LiS2 + L2S2 + L3S3)] . (B2) 

t=±i 

Here R is the normalization factor. 

If we take the zero temperature limit K ^ —00, the solution of the star-triangle relation 

is 



sinh 2Lj 



cosh2L, = 5L±fL^, ^2^2ziZ2^ 



k 



'2,ZjZf^ 



(B3) 



where 



and k) is a cychc permutation of (1,2,3). 

Now, we consider the transfer matrix T|f ^ on the honeycomb lattice whose matrix element 
T^^''(s,t) is the Boltzmann weight for a spin configuration shown in Fig. 12; 



T(f)(s,t)= Y: exp 



ri=±l 



X exp 



'M/2 

-^3^2m^2m+l 

m=l 



X exp 



M/2 



'M/2 

-^3^2m-l'"2n 

m=l 



X exp 



(-^1^2m-1^2m-l + -^2^2m^2m 

m=l 
■m/2 

(-^l'52m'"2m + -^2'52m-l'"2m-l 



m=l 



(T^ ■ Tbi ■ Tc ■ Tb2) is,t) 



(B4) 



The superscript fi (= 0, 1) denotes the boundary condition sm+ij = {—^Y^ij- Each of the 
four factors in the first two lines in Eq. ( [B^ defines the four matrices T^, T^i, Tc and Tb2, 
respectively. Their operator forms are 
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Um=i exp 



2m+l 



T*rT^ -X- T*rT 



X 

2m 



(B5) 



Ts2 = {^Y'^uflLW^ 

where erf is the Pauh spin operator at site i and e~^^» = tanhLj. Putting rows of Fig. 12 
in succession and applying the star-triangle transformation gives the M x N triangular 
lattice whose basis vectors are rotated by 90° from those of Fig. 1(a). Thus, to calculate 
the toroidal partition function of T = TAFIM, we instead carry out the calculation using 
Th'^^ The boundary condition along the vertical direction is Sj^jv+i = (— l)'^Sj_i. For each 
boundary condition, the partition function ^tapim is written as 



'^TAFIM 



Tr 



H 



AT 



(B6) 



where R is the spin reversal operator. 



Following the same procedure as in Ref . [p^ , we diagonalized the transfer matrix Eq, 



exactly. If we write the transfer matrix as Th'^ = (4zi2;2//i;^)'"'^ e then the hamiltonian 



Ti, after the usual Wigner- Jordan transformation can be written as 



1 



^ = 9 E Mk){2nk - 1) + e2(A;)(2mfc - 1)} 
^ k 



(B7) 



where and are the mutually commuting occupation number operators with eigenvalue 
Hk, rrik = 0, 1 and ei{k) and £2 (A;) are the quasi-particle excitation energy which is given by 



ei(fc) = sgn(/c) cosh ^t_|_ — cosh — icos ^ t_ 



e2{k) = sgn(A;) cosh ""^ it_|_ + cosh 



i cos ^ t 



where 



p 



1 

4 + zl 

2ZiZ2 



4 + 



^4 
^3 



Azlzl 



2z1 cos k 
Z1Z2 



(B9) 
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From now on, we set Z3 = 1 since all quantities are functions of zi/zs and -22/2^3 only. Zi/z^ 
in this section is equal to Xi of the text if £j = 0. (See Eq. (4).) The values of k are restricted 
to the set 



27T I 2Z - 1 + /i if Y.kiiT'k + rrik) = even 

/c = — X 

' 2Z + /i if 'Ek{nk + 'mk)= odd 

for boundary condition fi along the horizontal direction. 
The partition function ^tapiL is given by 

/ 9 X -Af/2 

^TAFIM ~ ^even ^ \ ^ ^oddj y^^J 



(BIO) 



(Bll) 



where 



even{odd) 



exp l-NY{ei{k)nk + e2{k)mk)'> 

'en(odd) \ k ) 



(B12) 



I^ovcn (Sodd) denoting the sums over the occupation number configurations {nk,mk = 0, 1} 
under the restriction J^ki'^k = even (odd), respectively, and the values of k are given 

in Eq. ( piO|) . We will say a state is in an even (odd) sector if Yliki^^k + ^fc) is even (odd). In 
Eq. ( [Bll| ) E^^ and are the ground state energy in the respective sectors and are given by 



E^ 



Et 



-9 L 1^1 

^ n=l 
^ n=l \ 



27r 



(2n-l+/i) 



M 

'2Tr{2n+iJ,j 
M 



+ 62 



27r 



+ 62 



M 

2ir{2n+fi) 
M 



{2n-l+fi) 



E 



ei 



27r 



M 

27r(2n+/i) 
M 



(2n-l+/i) 



(B13) 



We obtain the finite corrections to E^ and E^ from the Euler-Mclaurin formula. The results 
for periodic boundary condition (/i = 0) are 



E" 



Jo ^i{k)dk -jf (^j2 



a 



M_ 



J,'"^ e^{k)dk - '-^ {^^ - (1 - ay) ,«>i 



(B14) 



where 
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dk 



2iz! + zl)-l 



a 



Air 



kcM 
Air 



(B15) 



with [x] denoting the integer part of x. The ground state energies for each sector under 
anti-periodic boundary condition (/i = 1) are 



El = E° 



(B16) 



The quantity | In {Z3/2) — Jq"" ti{k)dk is the bulk free energy /oo per site. 

From the predictions of the conformal field theory, we know that the transfer matrix 
has gapless excitations with linear dispersion relation at the criticality. The quasi-particle 
excitation energies become zero aX k = ±kc where 

1 



cos kr 



2ZiZ2 



'z\ + zl) - {zl - zlf 



(B17) 



in the range \zi — < Z3 and \zi + 2:2! > ^3. So, we conclude that the system is critical in 
this range. This includes the result of Blote and Hilhorst p| who treated the case zi = Z2- 
The toroidal partition function .ZtapiIi under the general boundary condition (/i, v) is 



given as 



'TAFIM . , ^TAFIM/ 

M/N=fixed 



(B18) 



Especially, the toroidal partition function for periodic boundary condition in both directions 
is given as 



J(0,0) ^ Y S|«,^(i/i2-a2) -f4(l/12-(l/2-a)2)20 m^gN 

N/M=fixed 



Since N is large, the modes near k = ±kc whose energy scales like 1/M contribute factors of 
(9(1) in the sum. Therefore, for M, N 00 with N/M fixed, we may replace the dispersion 
relation by the linear one 



ei,2(fc) = fol^ =F kc\ - ia'^^k - k^) =F ia^ 



(B20) 
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where 



a 



da 
° " dk 



Zr, Z-l 



cosao = cos a (fee 



l-zl-zl 



2ZiZ2 

The restricted sums can be done conveniently using the transformation 

even 

E = {i + k^ {i - (-1)^^™'=} 



(B21) 



+ 



tr: 



(B22) 



where 



K= n E 

fc=^(2n-l)"fc = 

and similarly for tr^ and tr^. 



and 



tr" 



1 nt=0,l 



n E 

7 iirn ni.=0,l 



After a lengthy calculation with this linear dispersion relation Eq. 



^(0,0) 



{-|^i(^,g)p + |^2(^,g)P + |^3(^,g)P + l^4(^,g)P} 



, we obtain 

(B23) 



2|r/(g)P 

where r/(g) = g^/2^n^=i(l — g") is the Dedekind eta function, di are the Jacobi theta 
function |^ , q = e^'"*'^ with 



and finally 



^ = ^(eo + *«o) 
z = h vrar . 



Following the same procedures, we can also calculate ^tapim under the general boundary 
condition (/i, u). We present only the results. 

I;;|q2 

^(0,1) 



TAFIM 



2|r/(9)P 



2|^(g)P 
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(B24) 



The toroidal partition function >Z*'tafim can be rewritten as an infinite series in q. Using 
the series form of diS and rearranging the summands, we obtain 

^TAFiM — U(';^^|2 \ ^1 ^ 'i y \c>z,o) 

To compare with the triangular lattice shown in Fig. 1(a), we perform the modular 
transformation f— i>r = — 1/f. This is achieved by applying Poisson sum formula to both 
summation indices n and m in Eq. ( |B25| ). The resulting expression for the periodic boundary 
condition (yU = 0, z/ = 0) is 

7(0,0) _ ^ -2-wiam (:m/2+{n-aoN/2TT)f/2-{m/2-{n-aQN/2-K)Y/2 |'R9fi'l 



where q = e " . Note that 



with 



n-e, = cos-i 



and 



This is exactly the Coulombic partition function with the twisted boundary conditions and 
the coupling constant g = 2. 

Since Zj's are the activities of the diamonds shown in Fig. 1(b), one easily obtains p| 
from the bulk free energy that the mean domain wall densities of each type are given by 

(ni + n2)/N' = ao/iT 

{ni-n2)/M =kjTT (B27) 

with cto and are given by Eq. ( P21|) and ([B17| ), respectively. Therefore, one sees that 
the twisting angles are a^N and — 27ra in Eq. ( |B26| ) are related to the total domain wall 
densities as 

a„N = ^^<!^ . (B28) 
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and 

2.a = ^ (B29) 

where Qq is the number of domain walls per row and Qi is the difference of number of type 
1 domain walls and that of type 2 per column. Using this quantities, ^mriM is written as 



where /^m,n{9) = (jn/ ^ + ^n) /4. 
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Figure Captions 

Fig. 1. (a) Correspondence between the anisotropic couplings and the lattice directions in 
the TAFIM with nearest and next nearest neighbor interactions, (b) Labeling of three types 
of diamonds. Types 1 and 2 are considered as domain wall excitations. 
Fig. 2. A typical striped domain wall configuration derived from a ground state of the 
TAFIM on 4 X 6 lattice. The filled circles represent spin up states and empty circles represent 
spin down states. We use the periodic boundary condition along the horizontal direction and 
the anti-periodic boundary condition along the vertical direction. So, the resulting domain 
wall configuration has even Q and odd fi. 

Fig. 3. Four possibilities of next nearest neighbor bonds along the direction 1. 
Fig. 4. Deformation of Fig. 2 into a square lattice. 

Fig. 5. Five types of unit squares in the deformed lattice and assignment of vertex configu- 
rations. 

Fig. 6. Typical X2 dependence of q within the IC phase is shown for xi = X2- The curves 
are for A = —0.7, —0.3, 0, 0.3 and 0.7, respectively. 

Fig. 7. (a) Phase diagram in lna;2-A plane for a = 1. (0), (1) and (IC) denotes the C phase 
with g = 0, 1 and the IC phase, respectively, (b) Same as in (a) with a = 7. New C phase 
with q = 1/2 appears for A < —4. The Dotted line in the IC phase denotes the position 
where dislocations become irrelevant. See section |^ 

Fig. 8. A typical root distribution of the BAE in complex a plane. This figure shows the 
solutions of BAE for q = 1/2 and A = —7. 

Fig. 9. Typical X2 dependence of The domain wall densities for A = —4.0 and —5.0. The 
curves are for a = A + 0.1. This is obtained from the equation of state ( Eq. (p8D) where 
configurational free energy is taken from the numerical solution of the BAE with = 150. 
Fig. 10. The Gaussian coupling constant g calculated numerically with lattice size N = 150 
and A = 0.5, 0.0, —0.5, —2.0, —3.0, —4.0 and —5.0. For each curve, a is set to the value 
a = lAI +0.1. 
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Fig. 11. Effect of finite temperature in TAFIM is to excite dislocations in domain walls. 
Simultaneous creation of vortex (a) and anti-vortex (b) pairs as in (c) distroys the IC phase 
if 51 < 1 but is irrelevant if 51 > 1. 
Fig. 12. The honeycomb lattice transfer matrix. 
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